Instructions

Go through this R file in detail and understand the meaning of each script line. Answer the four questions (Q1 to Q4). Then organize the entire R file into a R Markdown file. Your R Markdown file should include each of the steps as below, including your answers to the four questions.

# set working directory
setwd("/Users/jonathanbernard/Desktop/UIUC/Spring 2024/GGIS 224/GGIS224/Labs/Lab10")

# load libraries
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(tmap)
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
library(leaflet)
# COUNTY BOUNDARY POLYGON DATA
# 2018 Illinois County Boundaries
# Source: Census, https://www.census.gov/geographies/mapping-files/time-series/geo/carto-boundary-file.html
counties.original <- st_read("Data/County Boundaries/cb_2018_us_county_500k.shp")
## Reading layer `cb_2018_us_county_500k' from data source 
##   `/Users/jonathanbernard/Desktop/UIUC/Spring 2024/GGIS 224/GGIS224/Labs/Lab10/Data/County Boundaries/cb_2018_us_county_500k.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 3233 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -179.1489 ymin: -14.5487 xmax: 179.7785 ymax: 71.36516
## Geodetic CRS:  NAD83
head(counties.original) # KEY = GEOID is 5-digit county ID 
## Simple feature collection with 6 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -89.18251 ymin: 36.91554 xmax: -83.45285 ymax: 38.52538
## Geodetic CRS:  NAD83
##   STATEFP COUNTYFP COUNTYNS       AFFGEOID GEOID    NAME LSAD      ALAND
## 1      21      007 00516850 0500000US21007 21007 Ballard   06  639387454
## 2      21      017 00516855 0500000US21017 21017 Bourbon   06  750439351
## 3      21      031 00516862 0500000US21031 21031  Butler   06 1103571974
## 4      21      065 00516879 0500000US21065 21065  Estill   06  655509930
## 5      21      069 00516881 0500000US21069 21069 Fleming   06  902727151
## 6      21      093 00516893 0500000US21093 21093  Hardin   06 1614569777
##     AWATER                       geometry
## 1 69473325 MULTIPOLYGON (((-89.18137 3...
## 2  4829777 MULTIPOLYGON (((-84.44266 3...
## 3 13943044 MULTIPOLYGON (((-86.94486 3...
## 4  6516335 MULTIPOLYGON (((-84.12662 3...
## 5  7182793 MULTIPOLYGON (((-83.98428 3...
## 6 17463238 MULTIPOLYGON (((-86.27756 3...
# Subset to just IL Data
counties.IL <- counties.original |> filter(STATEFP == 17) 
plot(counties.IL)

# ILLINOIS STATE BOUNDARY
il <- st_union(counties.IL) 
plot(il)

# COUNTY HEALTH ATTRIBUTE DATA
# 2023 Length of Life Estimates 
# Source: County Health Rankings, https://www.countyhealthrankings.org/explore-health-rankings/illinois/data-and-resources
# Technical Documentation: https://www.countyhealthrankings.org/sites/default/files/media/document/2023%20CHRR%20Technical%20Document.pdf
# Length of Life RANK = lower scores indicate BEST health, higher scores indicate WORSE health
lengthLife <- read.csv("Data/IL_LengthLife_CHR.csv")
head(lengthLife) #FIPS = County Key
##    FIPS    State    County LengthLife LLRank
## 1 17001 Illinois     Adams      -0.10     46
## 2 17003 Illinois Alexander       1.50    102
## 3 17005 Illinois      Bond      -0.13     40
## 4 17007 Illinois     Boone      -0.40     18
## 5 17009 Illinois     Brown      -0.80      4
## 6 17011 Illinois    Bureau      -0.16     38
# Merge Health Data to master county file
counties <- merge(counties.IL, lengthLife, by.x="GEOID", by.y="FIPS") #Essential a table join. 
head(counties) #Didn't work; inspect...CO_FIPS a bit strange!
## Simple feature collection with 6 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -91.51308 ymin: 36.9703 xmax: -88.70541 ymax: 42.49505
## Geodetic CRS:  NAD83
##   GEOID STATEFP COUNTYFP COUNTYNS       AFFGEOID      NAME LSAD      ALAND
## 1 17001      17      001 00424202 0500000US17001     Adams   06 2214804824
## 2 17003      17      003 00424203 0500000US17003 Alexander   06  609996948
## 3 17005      17      005 00424204 0500000US17005      Bond   06  985073312
## 4 17007      17      007 00424205 0500000US17007     Boone   06  727171389
## 5 17009      17      009 00424206 0500000US17009     Brown   06  791704469
## 6 17011      17      011 00424207 0500000US17011    Bureau   06 2250916543
##     AWATER    State    County LengthLife LLRank                       geometry
## 1 41767689 Illinois     Adams      -0.10     46 MULTIPOLYGON (((-91.51297 4...
## 2 44237171 Illinois Alexander       1.50    102 MULTIPOLYGON (((-89.51839 3...
## 3  6462629 Illinois      Bond      -0.13     40 MULTIPOLYGON (((-89.63926 3...
## 4  3361806 Illinois     Boone      -0.40     18 MULTIPOLYGON (((-88.94098 4...
## 5  4139310 Illinois     Brown      -0.80      4 MULTIPOLYGON (((-90.91703 3...
## 6 11606402 Illinois    Bureau      -0.16     38 MULTIPOLYGON (((-89.86235 4...
# Quick map for confirmation
# Length of Life RANK = lower scores indicate BETTER health, higher scores indicate WORSE health
tm_shape(il) + tm_borders(lwd = 2) +
  tm_shape(counties) + tm_polygons("LLRank", palette = "BuPu") +  
  tm_layout(frame = FALSE, legend.outside = TRUE)

# ADDING EMISSIONS POINT DATA
# 2020 National Emissions Inventory, Facility Data in IL
# Source: EPA, https://www.epa.gov/air-emissions-inventories/2020-national-emissions-inventory-nei-data

emissions.original <- read.csv("Data/NEI2020_FacilityIL.csv")
head(emissions.original)
##   pointID                              SITE..NAME EIS.Facility.ID    STATE
## 1       1                    ADM Animal Nutrition         3344211 Illinois
## 2       2                Prince Agri Products Inc         3345811 Illinois
## 3       3                     Prince Minerals Inc         7314711 Illinois
## 4       4 Quincy Regional Airport - Baldwin Field         3346911 Illinois
## 5       5                 Bunge North America Inc         3347011 Illinois
## 6       6                     Cairo Dry Kilns Inc         3362711 Illinois
##     State.County      Pollutant Pollutant.Type Emissions..Tons.
## 1     IL - Adams Lead Compounds        CAP/HAP          0.00001
## 2     IL - Adams Lead Compounds        CAP/HAP          0.00006
## 3     IL - Adams Lead Compounds        CAP/HAP          0.00256
## 4     IL - Adams Lead Compounds        CAP/HAP          0.03723
## 5 IL - Alexander Lead Compounds        CAP/HAP          0.00032
## 6 IL - Alexander Lead Compounds        CAP/HAP          0.00002
##                    Facility.Type               Street.Address
## 1                    Unspecified               1000 N 30th St
## 2                    Unspecified            4618 Gardner Expy
## 3                    Unspecified             401 N Prince Plz
## 4                        Airport             1645 Highway 104
## 5 Food Products Processing Plant                  203 34th St
## 6                    Unspecified 2 Miles N Of Cairo on Hwy 51
##                                     NAICS EPA.Region  FIPS Latitude Longitude
## 1         Other Animal Food Manufacturing          5 17001 39.94534 -91.36624
## 2         Other Animal Food Manufacturing          5 17001 39.87616 -91.39827
## 3 Synthetic Dye and Pigment Manufacturing          5 17001 39.89511 -91.40916
## 4                Other Airport Operations          5 17001 39.94093 -91.19425
## 5    Soybean and Other Oilseed Processing          5 17003 37.01605 -89.17777
## 6                            Home Centers          5 17003 37.04587 -89.18468
#Convert CSV to Spatial Data
emissions <- st_as_sf(emissions.original, 
                      coords = c("Longitude", "Latitude"),
                      crs = 4326)

# Check for type of pollutant in the dataset
unique(emissions$Pollutant)
## [1] "Lead Compounds"
# Using dots, adjust style, palette, and number of bins (n)
# By now, you are expected to be able to self-interpret what following statement does.
tm_shape(il) + tm_borders(lwd = 2) + #Draw state boundary
  tm_shape(counties) + tm_borders(lwd = 0.5) +  #overlay county boundaries
  tm_shape(emissions) +  #overlay point coordinates
  tm_dots("Emissions..Tons.", palette = "BrBG", #then apply graduated color based on attribute "Emissions..Tons."
          n = 10, style = "quantile",
          title="Lead Emissions (tons)") + 
  tm_layout(frame = FALSE, legend.outside = TRUE)

# Using bubbles, adjust size, color, style, and number of bins (n)
# Again, you are expected to be able to self-interpret what following statement does.
tm_shape(il) + tm_borders(lwd = 2) + #Draw state boundary
  tm_shape(counties) + tm_borders(lwd = 0.5) + #overlay county boundaries
  tm_shape(emissions) + #overlay point coordinates
  tm_bubbles("Emissions..Tons.", col = "tomato1") + #then apply graduated symbols based on attribute "Emissions..Tons."
  tm_layout(frame = FALSE, legend.outside = TRUE)

# Add the health data for exploration of associations with lead emissions.
tm_shape(il) + tm_borders(lwd = 2) +
  tm_shape(counties) + tm_polygons("LLRank", palette = "Greys") +  
  tm_shape(emissions) + 
  tm_bubbles(size = "Emissions..Tons.",  col = "Emissions..Tons.", 
             palette = "PuRd", style = "quantile") + 
  tm_layout(frame = FALSE, legend.outside = TRUE) 

tmap_mode("view")
## tmap mode set to interactive viewing
# Be sure to "Zoom" to explore map interactively
tm_shape(il) + tm_borders(lwd = 2) +
  tm_shape(counties) + tm_polygons("LLRank", alpha = 0.5, palette = "Greys") +  
  tm_shape(emissions) + 
  tm_bubbles(size = "Emissions..Tons.",  col = "Emissions..Tons.", 
             palette = "PuRd", style = "quantile") +  
  tm_basemap("Esri.WorldTopoMap")
## Legend for symbol sizes not available in view mode.
# Search for alternate basemaps -- 
# Use Leaflet Provider Demo at https://leaflet-extras.github.io/leaflet-providers/preview/
# & try out different Provider Names (in the tm_basemap argument)
tm_shape(il) + tm_borders(lwd = 2) +
  tm_shape(counties) + tm_polygons("LLRank", alpha = 0.5, palette = "Greys") +  
  tm_shape(emissions) + 
  tm_bubbles(size = "Emissions..Tons.",  col = "Emissions..Tons.", 
             palette = "PuRd", style = "quantile") +  
  tm_basemap("OpenStreetMap.HOT") 
## Legend for symbol sizes not available in view mode.
# Switch back to plot mode
tmap_mode("plot")
## tmap mode set to plotting
st_crs(emissions)
## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]
st_crs(counties)
## Coordinate Reference System:
##   User input: NAD83 
##   wkt:
## GEOGCRS["NAD83",
##     DATUM["North American Datum 1983",
##         ELLIPSOID["GRS 1980",6378137,298.257222101,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4269]]
emissions.3435 <- st_transform(emissions, 3435) #"NAD83 / Illinois East (ftUS)"
counties.3435 <- st_transform(counties, 3435)

st_crs(emissions.3435)
## Coordinate Reference System:
##   User input: EPSG:3435 
##   wkt:
## PROJCRS["NAD83 / Illinois East (ftUS)",
##     BASEGEOGCRS["NAD83",
##         DATUM["North American Datum 1983",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4269]],
##     CONVERSION["SPCS83 Illinois East zone (US Survey feet)",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",36.6666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-88.3333333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.999975,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",984250,
##             LENGTHUNIT["US survey foot",0.304800609601219],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["US survey foot",0.304800609601219],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["US survey foot",0.304800609601219]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["US survey foot",0.304800609601219]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["United States (USA) - Illinois - counties of Boone; Champaign; Clark; Clay; Coles; Cook; Crawford; Cumberland; De Kalb; De Witt; Douglas; Du Page; Edgar; Edwards; Effingham; Fayette; Ford; Franklin; Gallatin; Grundy; Hamilton; Hardin; Iroquois; Jasper; Jefferson; Johnson; Kane; Kankakee; Kendall; La Salle; Lake; Lawrence; Livingston; Macon; Marion; Massac; McHenry; McLean; Moultrie; Piatt; Pope; Richland; Saline; Shelby; Vermilion; Wabash; Wayne; White; Will; Williamson."],
##         BBOX[37.06,-89.28,42.5,-87.02]],
##     ID["EPSG",3435]]
st_crs(counties.3435) #102 counties
## Coordinate Reference System:
##   User input: EPSG:3435 
##   wkt:
## PROJCRS["NAD83 / Illinois East (ftUS)",
##     BASEGEOGCRS["NAD83",
##         DATUM["North American Datum 1983",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4269]],
##     CONVERSION["SPCS83 Illinois East zone (US Survey feet)",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",36.6666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-88.3333333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.999975,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",984250,
##             LENGTHUNIT["US survey foot",0.304800609601219],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["US survey foot",0.304800609601219],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["US survey foot",0.304800609601219]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["US survey foot",0.304800609601219]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["United States (USA) - Illinois - counties of Boone; Champaign; Clark; Clay; Coles; Cook; Crawford; Cumberland; De Kalb; De Witt; Douglas; Du Page; Edgar; Edwards; Effingham; Fayette; Ford; Franklin; Gallatin; Grundy; Hamilton; Hardin; Iroquois; Jasper; Jefferson; Johnson; Kane; Kankakee; Kendall; La Salle; Lake; Lawrence; Livingston; Macon; Marion; Massac; McHenry; McLean; Moultrie; Piatt; Pope; Richland; Saline; Shelby; Vermilion; Wabash; Wayne; White; Will; Williamson."],
##         BBOX[37.06,-89.28,42.5,-87.02]],
##     ID["EPSG",3435]]
nrow(emissions.3435) #850 emission facility points
## [1] 850
nrow(counties.3435) #102 counties
## [1] 102
# Count all of the facilities that intersect a county.
# we then use lengths() to find out how many items are present in a vector (i.e. county).
# Then, we store the facility counts for each of the counties into a new attribute called TotEmisFac.
counties.3435$TotEmisFac <- lengths(st_intersects(counties.3435, emissions.3435))

# Map the facility count per county as a choropleth map
tm_shape(il) + tm_borders(lwd = 2) +
  tm_shape(counties.3435) + tm_polygons(col = "TotEmisFac", n=10,
                                    style = "jenks", title = "Total Facilities") +
  tm_layout(frame = FALSE, legend.outside = TRUE)

# Calculate County Area
counties.3435$county_area <- st_area(counties.3435) 

Q1:

What is the unit of the calculated county area? Answer: the unit is square meters because the st_area() returns area values in square meters by default even though hte CRS used has units in US survey feet, the st_area() function still returns the area in square meters.

Q2:

Make a choropleth map showing the density of emission facilities (i.e., the count of emission points per square mile) by county.

counties.3435$density <- counties.3435$TotEmisFac / (as.numeric(counties.3435$county_area) / 27878400)

tm_shape(il) + tm_borders(lwd = 2) +  
  tm_shape(counties.3435) + tm_polygons(col = "density", n=10,                                    
                                    style = "jenks", title = "Facility Density (per sq mile)") +
  tm_layout(frame = FALSE, legend.outside = TRUE)

# Generate a point-in-polygon (PIP) dataset, joining county data to point data
pip <- st_join(emissions.3435, counties.3435, join = st_intersects)
head(pip)
## Simple feature collection with 6 features and 29 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 121147.9 ymin: 128297.8 xmax: 737702 ymax: 1208455
## Projected CRS: NAD83 / Illinois East (ftUS)
##   pointID                              SITE..NAME EIS.Facility.ID    STATE
## 1       1                    ADM Animal Nutrition         3344211 Illinois
## 2       2                Prince Agri Products Inc         3345811 Illinois
## 3       3                     Prince Minerals Inc         7314711 Illinois
## 4       4 Quincy Regional Airport - Baldwin Field         3346911 Illinois
## 5       5                 Bunge North America Inc         3347011 Illinois
## 6       6                     Cairo Dry Kilns Inc         3362711 Illinois
##     State.County      Pollutant Pollutant.Type Emissions..Tons.
## 1     IL - Adams Lead Compounds        CAP/HAP          0.00001
## 2     IL - Adams Lead Compounds        CAP/HAP          0.00006
## 3     IL - Adams Lead Compounds        CAP/HAP          0.00256
## 4     IL - Adams Lead Compounds        CAP/HAP          0.03723
## 5 IL - Alexander Lead Compounds        CAP/HAP          0.00032
## 6 IL - Alexander Lead Compounds        CAP/HAP          0.00002
##                    Facility.Type               Street.Address
## 1                    Unspecified               1000 N 30th St
## 2                    Unspecified            4618 Gardner Expy
## 3                    Unspecified             401 N Prince Plz
## 4                        Airport             1645 Highway 104
## 5 Food Products Processing Plant                  203 34th St
## 6                    Unspecified 2 Miles N Of Cairo on Hwy 51
##                                     NAICS EPA.Region  FIPS GEOID STATEFP
## 1         Other Animal Food Manufacturing          5 17001 17001      17
## 2         Other Animal Food Manufacturing          5 17001 17001      17
## 3 Synthetic Dye and Pigment Manufacturing          5 17001 17001      17
## 4                Other Airport Operations          5 17001 17001      17
## 5    Soybean and Other Oilseed Processing          5 17003 17003      17
## 6                            Home Centers          5 17003 17003      17
##   COUNTYFP COUNTYNS       AFFGEOID      NAME LSAD      ALAND   AWATER    State
## 1      001 00424202 0500000US17001     Adams   06 2214804824 41767689 Illinois
## 2      001 00424202 0500000US17001     Adams   06 2214804824 41767689 Illinois
## 3      001 00424202 0500000US17001     Adams   06 2214804824 41767689 Illinois
## 4      001 00424202 0500000US17001     Adams   06 2214804824 41767689 Illinois
## 5      003 00424203 0500000US17003 Alexander   06  609996948 44237171 Illinois
## 6      003 00424203 0500000US17003 Alexander   06  609996948 44237171 Illinois
##      County LengthLife LLRank TotEmisFac                    county_area
## 1     Adams       -0.1     46          4 24325392002 [US_survey_foot^2]
## 2     Adams       -0.1     46          4 24325392002 [US_survey_foot^2]
## 3     Adams       -0.1     46          4 24325392002 [US_survey_foot^2]
## 4     Adams       -0.1     46          4 24325392002 [US_survey_foot^2]
## 5 Alexander        1.5    102          3  7044304829 [US_survey_foot^2]
## 6 Alexander        1.5    102          3  7044304829 [US_survey_foot^2]
##       density                  geometry
## 1 0.004584247  POINT (133815.6 1208455)
## 2 0.004584247  POINT (123966.1 1183559)
## 3 0.004584247  POINT (121147.9 1190566)
## 4 0.004584247  POINT (181996.7 1205255)
## 5 0.011872740   POINT (737702 128297.8)
## 6 0.011872740 POINT (735782.5 139171.5)
# Using the PIP file, average all of emissions per county.
temp <- aggregate(pip$Emissions..Tons., by = list(pip$GEOID), mean)
head(temp)
##   Group.1            x
## 1   17001 9.965000e-03
## 2   17003 6.516667e-03
## 3   17005 1.269750e-02
## 4   17007 2.994800e-02
## 5   17009 3.390000e-03
## 6   17011 2.666667e-05
# Rename the data in temp file
names(temp) <- c("GEOID", "EmissionAve")

# Join back to master county dataset
counties.3435 <- merge(counties.3435, temp, by = "GEOID", all.x = TRUE)
head(counties.3435)
## Simple feature collection with 6 features and 17 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 95633.61 ymin: 111619.1 xmax: 883631.1 ymax: 2123580
## Projected CRS: NAD83 / Illinois East (ftUS)
##   GEOID STATEFP COUNTYFP COUNTYNS       AFFGEOID      NAME LSAD      ALAND
## 1 17001      17      001 00424202 0500000US17001     Adams   06 2214804824
## 2 17003      17      003 00424203 0500000US17003 Alexander   06  609996948
## 3 17005      17      005 00424204 0500000US17005      Bond   06  985073312
## 4 17007      17      007 00424205 0500000US17007     Boone   06  727171389
## 5 17009      17      009 00424206 0500000US17009     Brown   06  791704469
## 6 17011      17      011 00424207 0500000US17011    Bureau   06 2250916543
##     AWATER    State    County LengthLife LLRank TotEmisFac
## 1 41767689 Illinois     Adams      -0.10     46          4
## 2 44237171 Illinois Alexander       1.50    102          3
## 3  6462629 Illinois      Bond      -0.13     40          4
## 4  3361806 Illinois     Boone      -0.40     18          5
## 5  4139310 Illinois     Brown      -0.80      4          2
## 6 11606402 Illinois    Bureau      -0.16     38          3
##                      county_area     density  EmissionAve
## 1 24325392002 [US_survey_foot^2] 0.004584247 9.965000e-03
## 2  7044304829 [US_survey_foot^2] 0.011872740 6.516667e-03
## 3 10673939103 [US_survey_foot^2] 0.010447277 1.269750e-02
## 4  7864025133 [US_survey_foot^2] 0.017725274 2.994800e-02
## 5  8576601423 [US_survey_foot^2] 0.006501037 3.390000e-03
## 6 24360918399 [US_survey_foot^2] 0.003433171 2.666667e-05
##                         geometry
## 1 MULTIPOLYGON (((95734.59 12...
## 2 MULTIPOLYGON (((639491.9 22...
## 3 MULTIPOLYGON (((613093.3 85...
## 4 MULTIPOLYGON (((819702.8 20...
## 5 MULTIPOLYGON (((259336 1188...
## 6 MULTIPOLYGON (((565941.6 17...
# Map EmissionAve
tm_shape(il) + tm_borders(lwd = 2) +
  tm_shape(counties.3435) + tm_polygons(col = "EmissionAve", n=10,
                                    style = "jenks", title = "Emission Ave") +
  tm_layout(frame = FALSE, legend.outside = TRUE) 

# Generate 5-mile buffers. Note: 1 mile = 5280 ft (using EPSG 3435)
emission_buffers <- st_buffer(emissions.3435, 5280*5)
glimpse(emission_buffers) #Take a look at the transposed data: check https://dplyr.tidyverse.org/reference/glimpse.html
## Rows: 850
## Columns: 14
## $ pointID          <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16…
## $ SITE..NAME       <chr> "ADM Animal Nutrition", "Prince Agri Products Inc", "…
## $ EIS.Facility.ID  <int> 3344211, 3345811, 7314711, 3346911, 3347011, 3362711,…
## $ STATE            <chr> "Illinois", "Illinois", "Illinois", "Illinois", "Illi…
## $ State.County     <chr> "IL - Adams", "IL - Adams", "IL - Adams", "IL - Adams…
## $ Pollutant        <chr> "Lead Compounds", "Lead Compounds", "Lead Compounds",…
## $ Pollutant.Type   <chr> "CAP/HAP", "CAP/HAP", "CAP/HAP", "CAP/HAP", "CAP/HAP"…
## $ Emissions..Tons. <dbl> 0.00001, 0.00006, 0.00256, 0.03723, 0.00032, 0.00002,…
## $ Facility.Type    <chr> "Unspecified", "Unspecified", "Unspecified", "Airport…
## $ Street.Address   <chr> "1000 N 30th St", "4618 Gardner Expy", "401 N Prince …
## $ NAICS            <chr> "Other Animal Food Manufacturing", "Other Animal Food…
## $ EPA.Region       <int> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,…
## $ FIPS             <int> 17001, 17001, 17001, 17001, 17003, 17003, 17003, 1700…
## $ geometry         <POLYGON [US_survey_foot]> POLYGON ((160215.6 1208455,..., …
# Map the buffers
tm_shape(il) + tm_borders(lwd = 2) +
  tm_shape(counties) + tm_polygons(alpha = 0.5) +  
  tm_shape(emission_buffers) + tm_borders(col = "tomato1") +
  tm_shape(emissions.3435) + tm_dots(col = "black")

########################
# CALCULATE BUFFER COUNT BY COUNTY
########################

# First we identify how many times buffers overlap counties.
# Then we use lengths() to find out how many items are present in each county. 
counties.3435$bufferCt <- lengths(st_intersects(counties.3435, emission_buffers))
head(counties.3435)
## Simple feature collection with 6 features and 18 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 95633.61 ymin: 111619.1 xmax: 883631.1 ymax: 2123580
## Projected CRS: NAD83 / Illinois East (ftUS)
##   GEOID STATEFP COUNTYFP COUNTYNS       AFFGEOID      NAME LSAD      ALAND
## 1 17001      17      001 00424202 0500000US17001     Adams   06 2214804824
## 2 17003      17      003 00424203 0500000US17003 Alexander   06  609996948
## 3 17005      17      005 00424204 0500000US17005      Bond   06  985073312
## 4 17007      17      007 00424205 0500000US17007     Boone   06  727171389
## 5 17009      17      009 00424206 0500000US17009     Brown   06  791704469
## 6 17011      17      011 00424207 0500000US17011    Bureau   06 2250916543
##     AWATER    State    County LengthLife LLRank TotEmisFac
## 1 41767689 Illinois     Adams      -0.10     46          4
## 2 44237171 Illinois Alexander       1.50    102          3
## 3  6462629 Illinois      Bond      -0.13     40          4
## 4  3361806 Illinois     Boone      -0.40     18          5
## 5  4139310 Illinois     Brown      -0.80      4          2
## 6 11606402 Illinois    Bureau      -0.16     38          3
##                      county_area     density  EmissionAve
## 1 24325392002 [US_survey_foot^2] 0.004584247 9.965000e-03
## 2  7044304829 [US_survey_foot^2] 0.011872740 6.516667e-03
## 3 10673939103 [US_survey_foot^2] 0.010447277 1.269750e-02
## 4  7864025133 [US_survey_foot^2] 0.017725274 2.994800e-02
## 5  8576601423 [US_survey_foot^2] 0.006501037 3.390000e-03
## 6 24360918399 [US_survey_foot^2] 0.003433171 2.666667e-05
##                         geometry bufferCt
## 1 MULTIPOLYGON (((95734.59 12...        4
## 2 MULTIPOLYGON (((639491.9 22...        3
## 3 MULTIPOLYGON (((613093.3 85...        7
## 4 MULTIPOLYGON (((819702.8 20...       19
## 5 MULTIPOLYGON (((259336 1188...        3
## 6 MULTIPOLYGON (((565941.6 17...       19
tm_shape(il) + tm_borders(lwd = 2) +
  tm_shape(counties.3435) + tm_polygons("bufferCt", n=10,
                    style = "jenks", title = "CountBuffer") + 
  tm_layout(frame = FALSE, legend.outside = TRUE)

# Q3: Compare the buffer count by county vs the point count by county. What difference do you notice? Answer: The buffer count by country is always much higher than the point count by county. Although the an emission point is only counted once per exact location, the 5-mile buffer around each point can interserct multiple counties. Hence, leading to a much higher count of affected counties when using buffer zones compared to just using the point locations.

Q4: Further quantify the effects of emission diffusion/dispersion

Integrate what you have learned thus far to answer the following two questions. There are several ways to solve each of the questions.

Q4.1. Calculate the number of counties affected by each of the emission buffer zones.

Write your code chunk below and map the result as a GRADUATED SYMBOLOGY MAP. Be creative in the mapping design.

emission_buffers$affectedCounties <- lengths(st_intersects(emission_buffers, counties.3435))

tm_shape(il) + tm_borders(lwd = 2) + 
  tm_shape(emission_buffers) +
  tm_polygons(col = "affectedCounties",             
              palette = "YlOrRd", style = "quantile", title = "Affected Counties",
              legend.is.portrait = TRUE) +
  tm_layout(frame = FALSE, legend.outside = TRUE)

Q4.2. Calculate the percentage of area in each county affected by emission diffusion.

Write your code chunk below and map the result as a choropleth MAP.

intersection <- st_intersection(counties.3435, emission_buffers)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
pctAffected <- aggregate(as.numeric(st_area(intersection)), by = list(intersection$GEOID), FUN = sum)
names(pctAffected) <- c("GEOID", "affectedArea")

counties.3435 <- merge(counties.3435, pctAffected, by = "GEOID", all.x = TRUE)
counties.3435$pctAffected <- counties.3435$affectedArea / as.numeric(counties.3435$county_area) * 100

tm_shape(il) + tm_borders(lwd = 2) +
  tm_shape(counties.3435) + tm_polygons(col = "pctAffected", n = 10,
                                        style = "jenks", title = "% Area Affected",
                                        legend.is.portrait = TRUE) + 
  tm_layout(frame = FALSE, legend.outside = TRUE)